Inertial Frame Independent Forcing for Discrete Ve- 
locity Boltzmann Equation: Implications for Filtered 
Turbulence Simulation 

Kannan N. Premnath 1 '*, and Sanjoy Banerjee 2 

1 Department of Mechanical Engineering, University of Wyoming, Laramie, WY 
82071, U.S.A. 

2 Department of Chemical Engineering, City College of New York, City University of 
New York, New York, NY 10031, U.S.A. 



Abstract. We present a systematic derivation of a model based on the central moment 
lattice Boltzmann equation that rigorously maintains Galilean invariance of forces to 
simulate inertial frame independent flow fields. In this regard, the central moments, 
<-j^*! i-e- moments shifted by the local fluid velocity, of the discrete source terms of the lattice 

Boltzmann equation are obtained by matching those of the continuous full Boltzmann 
equation of various orders. This results in an exact hierarchical identity between the 
central moments of the source terms of a given order and the components of the central 
moments of the distribution functions and sources of lower orders. The corresponding 
source terms in velocity space are then obtained from an exact inverse transformation 
due to a suitable choice of orthogonal basis for moments. Furthermore, such a cen- 
^ I tral moment based kinetic model is further extended by incorporating reduced com- 

p> ' pressibility effects to represent incompressible flow. Moreover, the description and 

simulation of fluid turbulence for full or any subset of scales or their averaged be- 
havior should remain independent of any inertial frame of reference. Thus, based on 
the above formulation, a new approach in lattice Boltzmann framework to incorporate 
turbulence models for simulation of Galilean invariant statistical averaged or filtered 

turbulent fluid motion is discussed. 
PACS: 05.20.Dd,47.27.-i,47.27.E- 
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1 Introduction 

Minimal kinetic models for the Boltzmann equation, i.e. lattice Boltzmann equation for- 
mulations, are evolving towards as alternative physically-inspired computational tech- 
niques for various fluid mechanics and other problems. Originally developed as an im- 
proved variant of the lattice gas automata HI to eliminate statistical noise 0, the lattice 
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Boltzmann method (LBM) has undergone a series of major refinements, in terms of its 
underlying physical models as well as numerical solution schemes for various applica- 
tions over the last two decades |3H6]|. In particular, its rigorous connection to the kinetic 
theory IZrlSl has resulted in a number of recent developments, including models that 
are more physically consistent for multiphase l ilOllTTII and multicomponent flows |12|, 
models for non-equilibrium phenomena beyond the Navier-Stokes-Fourier representa- 
tion [ 13 1 and an asymptotic analysis approach to establish consistency of the LBM from a 
numerical point of view |14| . 

The stream-and-collide procedure of the LBM can be considered as a dramatically 
simplified discrete representation of the continuous Boltzmann equation. Here, the stream- 
ing step represents the advection of the distribution of particle populations along discrete 
directions, which are designed from symmetry considerations, between successive colli- 
sions. Much of the physical effects being modeled are represented in terms of the collision 
step, which also significantly influences the numerical stability of the LBM. Most of the 
major developments until recently were associated with the single-relaxation-time (SRT) 
model 1 15, 16 1 based on the BGK approximation fTTfl, and enjoys its popularity owing, 
mainly, to its simplicity. However, it is prone to numerical instability. Moreover, it is 
inadequate in its representation of certain physical aspects, such as independently ad- 
justable transport properties of thermal transport and viscoelastic phenomena. 

These limitations have been significantly addressed in the multiple-relaxation-time 
(MRT) collision model IT8l . This, in a sense, represents a simplified form of the relaxation 
LBM proposed earlier ||19U20|) , with an important characteristic difference in that the col- 
lision process is carried out in moment space |2TH instead of in the usual velocity space. 
By separating the relaxation time scales of different moments, obtained by using a linear 
Fourier stability analysis, its numerical stability can be significantly improved Il22ll23l . 
Furthermore, it has resulted in significant advantages over the SRT-LBM for computa- 
tion of various classes of fluid flow problems, including multiphase systems [24-26], tur- 
bulent flows [27,28] and magnetohydrodynamics 1129 II . It may be noted that recently a 
different form of MRT model based on the orthogonal Hermite polynomial projections 
of the distribution functions, which is independent of any underlying lattice structure, 
allowing representation of higher order non-equilibrium effects has been proposed [30]. 

The stabilization of the LBM using a single relaxation time has been addressed from a 
different perspective by enforcing the H-theorem locally in the collision step (3T14341 . By 
using the attractors of the distribution function based on the minimization of a Lyapunov- 
type functional, non-linear stability of the LBM is achieved in this Entropic LBM. This ap- 
proach has recently been significantly extended to incorporate multiple relaxation times 
with efficient implementation strategies 155113611 . Furthermore, systematic procedures for 
different types of higher-order LBM have been developed |37l - [39| . An important element 
is the construction of higher-order lattices based on symmetry considerations which have 
been analyzed using group theory [40,41]. Further progress, from a numerical aspect, is 
that based on the consistency analysis [14| and a notion of structural stability |[42ll43l 
(shown related to the Onsager-like relation in non-equilibrium thermodynamics [44]), 
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convergence of the LBM to the Navier-Stokes equations has rigorously been shown [45 J. 

On the other hand, it is important to clearly understand in what sense the lattice Boltz- 
mann equation (LBE), which is generally considered as a mesoscopic approach, inherits 
or maintains the various physical invariance properties of the continuous full Boltzmann 
equation (which it represents as a much simplified model) and the Navier-Stokes equa- 
tions (which it represents numerically). Careful considerations of these aspects play an 
important role in ensuring the general applicability of the approach for various, espe- 
cially challenging, problems. In this regard, and to put the present work in perspective, 
it should be noted that the continuum mechanics description as well as the microscopic 
statistical (continuous Boltzmann) description of fluid motion generally satisfy a larger 
invariance group, with inertial frame invariance being just an important special case. 
The most general form among these is the so-called the principal of material frame indif- 
ference, also known as the objectivity principle [46 J. According to this, the constitutive 
equations should have the same forms in all frames of reference, whether inertial or not. 
While this is considered as an important axiom based on which the continuum mechanics 
is formulated 06), its role from continuous kinetic theory point of view was the subject 
of considerable analysis for sometime |[47l - l52| . 

The following are the main outcomes of these studies: the continuous full Boltzmann 
equation (i) is material frame indifferent in a strong approximate sense, when there is a large 
scale separation between the collision times and the macroscopic flow times B9ll5Tll52| 
(thus providing a strong support to the axiomatic principle generally used in the contin- 
uum description), (ii) satisfies both the inertial frame or Galilean invariance as well as 
the extended Galilean invariance (i.e. invariance under arbitrary translational accelera- 
tions of the reference frame) exactly [52J. Furthermore, it was shown that the standard 
procedures (e.g. Chapman-Enskog expansion, Maxwellian iteration) lead to frame de- 
pendent higher order contributions for the constitutive equations in non-inertial frames 
in the continuum limit, while the continuous kinetic theory itself can be frame indepen- 
dent Il49ll . Careful considerations of these principles could guide in the development of 
more generally applicable models and numerical schemes for complex problems. For ex- 
ample, material frame indifference (point (i)) is generally used as an important constraint 
for the constitutive equations for complex fluids (e.g. beyond Newtonian constitutive 
description such as polymeric fluids) and in the development of turbulence models in 
continuum mechanics. As mentioned above, this property is satisfied in a strong ap- 
proximate sense by the continuous Boltzmann equation, but not necessarily by the tools 
that relate the microscopic and macroscopic descriptions(point iii). This aspects are per- 
tinent in the construction of complex models from the continuous Boltzmann equation 
(e.g. [53, 54 1). In this work, however, we limit our discussion to an association of the 
properties mentioned in a part of the point (ii) for the LBE, i.e. for the exact invariance 
group - the Galilean or inertial frame invariance. 

In this regard, as a discrete approximation to the continuous full Boltzmann equation, 
the development of the LBE consists of simplifications at different levels. Thus, its vari- 
ous elements should be analyzed carefully to ascertain and quantify as to how well it sat- 
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isfies Galilean invariance. First, in contrast to continuous kinetic theory, due to the choice 
of finite lattice velocity sets and associated symmetries, it introduces linear dependencies 
of higher order moments with those of lower order moments that are supported by the 
lattice set |40| . Such degeneracies can in turn lead to negative dependence of viscosity on 
fluid velocity. It generally causes the Galilean invariance to be broken by the presence of 
terms that are cubic in velocity for the standard lattice configurations (with symmetries of 
square in 2D and cube in 3D) and also leads to numerical instability, especially at higher 
Mach numbers. This issue can be alleviated by the use of extended lattice velocity sets, 
which then relegates the degeneracies among moments to even higher orders. Second, 
the collision step including the forcing terms of the LBE should be carefully constructed 
in such a way that they recover correct physics which is inertial frame independent, i.e. 
the Navier-Stokes equations. Here, the use of independent set of central moments for a 
chosen lattice provides a natural approach to maintain Galilean invariance that can be 
constructed by invoking elements directly from kinetic theory. This is the main goal of 
the present work (see below). A rational means to more efficiently account for both the 
above aspects is discussed in the last section of this paper. And, third, the streaming step 
of the LBE is generally constructed as a discrete Lagrangian process. In the standard im- 
plementation, this couples the particle velocity and configurations spaces, which in turn, 
constrains the numerical accuracy of the LBE in the representation of the Navier-Stokes 
equations. As a result, the Galilean invariance of the LBE is limited by its overall nu- 
merical accuracy. However, it is known that such coupling between physical and lattice 
symmetries is not necessary in the discretization of the streaming operator. In fact, it can 
be discretized using classical schemes such as finite-difference or finite-element methods 
that alleviate this issue (55l - l57ll . Specifically, exploiting higher order discretization and 
time integration schemes (e.g. 115810 for the streaming operator could further improve the 
order of accuracy (and hence the Galilean invariance) of the LBE. Furthermore, the use of 
implicit schemes could enhance the computational efficiency in this regard. 

Focusing on the second aspect mentioned above, a different type of collision opera- 
tor and forcing can be devised that can maintain Galilean invariance for a chosen lattice 
velocity set and a discretization scheme for the streaming step. Specifically, central mo- 
ments are relaxed in a moving frame of reference during collision step [59|, originally 
proposed to improve numerical stability, but emphasized here for its better physical co- 
herence. The use of central moments, which are obtained by shifting the particle velocity 
with the local fluid velocity [60J, rigorously enforces Galilean invariance. In particular, 
while other previous approaches are generally Galilean invariant for up to second-order 
moments, the central moment based approach provides a higher order frame invariance 
as permitted by the discrete lattice velocity set. This approach was examined based on 
the concept of generalized local equilibrium |6Tfl . In addition, to further improve phys- 
ical coherence, the attractors for the higher order central moments were constructed as 
products of the lower order central moments, leading to the factorized central moment 
method \62\. Recently, a new approach to incorporate source terms using central mo- 
ments in the LBM that are Galilean invariant by construction, which are important for 
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computation of various physical problems, was developed 163 j. The consistency of this 
technique to the Navier-Stokes equations was shown by means of the Chapman-Enskog 
analysis [64| and its numerical accuracy was established. Furthermore, the method was 
also extended in three-dimensions for various lattice velocity sets and validated for a 
class of canonical problems |65|. As clarified in [63,65], numerical stability of the central 
moment approach can be enhanced, when it is executed in a multiple relaxation time 
formulation, similar to the standard or raw moment based approaches. Interestingly, it 
has been shown recently that when some classical schemes for flow simulation are made 
to satisfy Galilean invariance more rigorously, they led to more robust implementations 

(e.g. ESUS3). 

Turbulence remains as among the most challenging classes of flows for which con- 
siderable effort has been focused on the development of theory and applications using 
the LBM. Since its roots can be traced to kinetic theory, the LBM has been analyzed 
for the development of turbulence models from a fundamental point of view II68H70II . 
It has been employed for computation of Reynolds-averaged description of turbulent 
flows H7TII72II . Furthermore, it has found applications for large eddy simulation (LES) 
using LBM formulations with SRT [73j, and MRT [74| with multiblock approach for ef- 
ficient implementation 1I27II . Recently, dynamic subgrid scale (SGS) models for LES were 
incorporated into the LBM framework that resulted in reduced empiricism for descrip- 
tion at such scales |28]| . Moreover, an improved inertial-range consistent SGS model was 
also proposed II75II . A theoretical formulation for a SGS model based on an approximate 
deconvolution method [78] that does not rely on the common eddy-viscosity concept for 
application with the LBM was also devised recently II76II . Lastly, the closure modeling 
issues of kinetic and continuum turbulence effects were reconciled in a unified statisti- 
cal/ filtered description using a modified kinetic equation 11771 . Effectively, this allows the 
use of macroscopic turbulence models involving divergence of the Reynolds stress in the 
forcing term of the kinetic equation. 

An important physical consideration for any description of turbulent flow is that it 
should be invariant for all inertial frames of reference. In other words, for general appli- 
cability, representation of turbulence for all or any subset of its scales should be Galilean 
invariant. Thus, in particular, all SGS models, and associated numerical schemes for 
turbulence computation, should be frame invariant. An insightful analysis of various 
turbulence models was carried out from this viewpoint in [79]. A method to achieve 
Galilean invariance by means of certain redefinition of turbulent stresses was discussed 
in |80| . A recent review on this subject is reported in [81J. Furthermore, it should be noted 
that concepts based on central moments have played an important role in the develop- 
ment of theoretical foundation of turbulence physics - such as for statistical turbulence 
theory [82 j and turbulence modeling [83 J. 

Thus, in this paper, we develop a lattice Boltzmann equation based on central mo- 
ments for Galilean invariant representation of turbulent flows. Specifically, it allows 
frame-independent incorporation of general models for turbulent Reynolds stresses in 
a statistical/ filter averaged formulation using LBM for turbulence simulation. Further- 
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more, in a general setting, it maintains the forces and stresses to be independent of any 
inertial frame of reference and could also improve numerical stability in computations. 
In 1 631, we developed a forcing scheme based on a particular ansatz involving the local 
Maxwell distribution. Here, we develop a general forcing based on central moments by 
a direct examination of the continuous full Boltzmann equation itself, which unlike 1 63 1 
could also self-consistently account for non-equilibrium effects in higher order terms. 
In this regard, the central moments of the resulting source terms of the continuous and 
discrete counterparts are matched successively at different orders leading to a cascaded 
structure. In essence, this approach can be considered as a Galilean invariant minimal 
discrete model for the full Boltzmann equation including forcing terms. The attr actors 
for higher order central moments in the collision step of this computational model is 
based on the factorization in terms of those at lower orders by including such general 
forcing terms. In addition, we further develop this approach with reduced compressibil- 
ity effects for improved representation of turbulent flow physics in the incompressible 
limit. The forcing formulation developed here for incorporating turbulence models in a 
statistical/ filtered formulation can be extended to other problems, such as, for example, 
Galilean invariant representation of forces or stresses in complex fluids. 

The paper is organized as follows. Section [2] briefly discusses the choice of the mo- 
ment basis employed in this paper and Sec. [3] the continuous Boltzmann equation. In 
Sees. |H and |5l continuous forms of the central moments for the distribution functions 
and its local equilibrium, and sources due to force fields, respectively, are introduced. 
The LBE based on central moments with the general forcing terms is presented in Sec. 
Various discrete central moments are presented in Sec. that also specifies a matching 
principle to maintain Galilean invariance and the relationships among such moments are 
provided in Sec. [8j Section [9] describes various discrete raw moments and the derivation 
of the source terms in terms of the discrete particle velocity space. In Sec.[lOlwe present 
the construction of the collision operator of the central moment based LBM. The compu- 
tational procedure of this approach is provided in Sec. [TTJ The derivation is extended by 
considering reduced compressibility effects in SecHH Furthermore, Sec.[l3]discusses the 
use of attractors of the higher order central moments based on the concept of their factor- 
ization in term of those at lower orders. A natural consequence of this overall approach 
is that turbulence models can be represented for Galilean invariant filtered turbulence 
simulation using the LBM, which is described in Sec. HH Finally, the summary and con- 
clusions of this work are discussed in Sec.[l5l 



2 Selection of Moment Basis 

An important element in the development of the central moment based LBM is the spec- 
ification of a suitable basis for moments. In this work, to elucidate our approach, the 
two-dimensional, nine velocity (D2Q9) lattice model (see Fig. [1} is considered, for which 
the moment basis used in |63l| is adopted. It should, however, be noted that the procedure 
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described henceforth is of general nature, and can be extended for other lattice models 
and in three dimensions. The particle velocity for this lattice model e a is given by 

C (0,0) a = 

e>=\ (±1,0),(0,±1) « = l / - / 4 (2.1) 
[ (±1,±1) a = 5, — ,8 

For convenience, we employ Dirac's bra-ket notion to represent the basis vectors, and 
Greek and Latin subscripts for particle velocity directions and Cartesian coordinate di- 
rections, respectively. Noting that moments in the LBM are discrete integral properties 
of the distribution function f a , i.e. Tl^o e Tx e "yfa, where m and n are integers, we begin 
with the following nine non-orthogonal independent basis vectors obtained by combin- 
ing monomials e™ x e" in an ascending order. That is, 



-11^1°) 


= (1,1,1,1,1,1,1,1,1) T , 




= (0,1,0,-1,0,1,-1,-1,1) T , 




= (0,0,1,0,-1,1,1,-1,-1) T , 


\e 2 +e 2 ) 


= (0,1,1,1,1,2,2,2,2) T , 


\e 2 -e 2 ) 

\^ax ay/ 


= (0,1, -1,1, -1,0,0,0,0) T , 


l^cex^ixy) 


= (0,0,0,0,0,1, -1,1, -1) T , 


\ e ax e °ty) 


= (0,0,0,0,0,1,1,-1,-1) T , 


\ e ax^\y) 


= (0,0,0,0,0,1,-1,-1,1) T , 


\e 2 e 2 ) 


= (0,0,0,0,0,1,1,1,1) T , 
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where the superscript 'T represents the transpose operator. 

For an efficient implementation, the above non-orthogonal basis set is transformed 
into an equivalent orthogonal set through the Gram-Schmidt procedure in the increas- 
ing order of the monomials of the products of the Cartesian components of the particle 
velocities 163: 



l*o> = 


= \p), 


l*i> = 


— | &KX ) / 


I-K2) = 


— £(XV ) f 
1 1 


l*3> = 


- ^I/J 2 4-f 2 \ 


\K A ) -- 


= le 2 -e 2 ) 


\K 5 ) - 




\Ke) ~- 


- 3 If 2 p S -1-2 If ) 


\K 7 ) 


= -3\e n ve 2 )+2\e„r) 


\Ks) - 


= 9|4e 2 y )-6|e 2 x +^ y )+4|p). 


This can be written explicitly in 


term of a matrix given by 




- 1 


-4 4 ■ 




1 1 


0-1 1 2-2 




1 


1 -1 -1 2 -2 




1 -1 


0-1 1 0-2-2 


K = 


1 


-1 -1 -1 0-2 -2 




1 1 


1 2 1-1-1 1 




1 -1 


1 2 0-1-1 1 1 




1 -1 


-12 1111 




. 1 1 


-1 2 0-1 1-1 1 . 


where we have used 






JC = 


[|Ko),|Ki),|K 2 >, \K 3 ), |X 4 ), |X 5 ), \K 6 ),\K 7 ),\K 8 )] 



(2.3) 



(2.4) 



(2.5) 



3 Continuous Boltzmann Equation 

We consider the two-dimensional (2D) continuous Boltzmann equation, for which we 
aim to develop a Galilean invariant discrete model using the above basis vectors. It rep- 
resents the evolution of the continuous density distribution function f = f(x,y,^ x ,^y) in 

continuous phase space (x,y,^ x ,^y) subjected to a local force field F = (F x ,Fy), whose 
origin could be internal or external to the system. By definition, the averaged effects of /, 
weighted by various powers of the continuous particle velocity (£ x £y)i i.e. its moments, 
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are considered to characterize the various physical processes inherent in the motion of 
athermal fluids. In particular, the evolution of the slow hydrodynamical processes are 
described by the local macroscopic fluid density p and fluid velocity u = (u x ,Uy). The 
continuous Boltzmann equation may be written as 



dt 
where 



9/ + f • — • ^ t f=n(f,f), (3.i) 



CO /'CO 



P = fd&dty, (3.2) 



OO J —CO 
OO f oo 



pH = fgdgA. (3.3) 



y 

— oo 7 —oo 



Here, Q (/,/) is the collision term, which represents the cumulative effect of binary col- 
lision of particles. The force fields modify the distribution function exactly by the term 

— p ~^~?f' wm ch is obtained by moving the last term on the left hand side of Eq. (13.11) 
to its right to serve as a source term. It was shown by Grad (1949) |2T| that that solu- 
tion of Eq. d3.1|) can be approximated by the evolution equations for a hierarchical set of 
moments. Here, we seek to obtain a dramatically discretized version of this continuous 
Boltzmann equation by means of a hierarchy of central moments, focusing, in particular, 
on the forcing term, to obtain Galilean invariant representation of the dynamics of fluid 
motion. 



4 Continuous Central Moments: Distribution Function and its 
Local Attractor 



We now consider the integral properties of the distribution function / given in terms of 
its central moments, i.e. those shifted by the macroscopic fluid velocity. In particular, we 
define continuous central moment of / of order (tn+n) as 



f(Z x -u x ) m (£ y -u y ) n dZ x dZ r 



(4.1) 



— CO J — oo 



Here, and in the rest of this paper, the use of "hat" over a symbol represents quantities 
in the space of moments. The distribution function for an athermal fluid has a local 
equilibrium state in the continuous particle velocity space (£ x >£y), which is given by the 
Maxwellian as |84l 



2d 



(4.2) 
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where Cg = 1/3. Analogously, we can define the corresponding central moment of the 
Maxwell distribution of order (m+n) as 

/CO /-CO 
/ f M {Z X -U X T{Z y -U y YdS X dZ y . (4.3) 
-CO J —CO 

By virtue of the fact that f M being an even function, Tl^y* ^0 when m and n are even and 
n^ v „ =0 when morn odd. Here and henceforth, the subscripts x m y n mean xxx ■ ■ ■ m-times 
and yyy ■ ■ ■ n -times. Evaluation of the central moments of the Maxwellian, to different or- 
ders of increasing powers, yields 

n<T = P , 
= o, 



M 
y 

M _ „2, 



nf = o, 



±± xx u sr' 

n$ = <* P , (4.4) 

= o, 

= o 

n M = o 

1 xyy u / 
±1 xxyy L sr - 

5 Continuous Central Moments: Forcing 

In the presence of a force field F, in view of Eq. d3.ll) and as discussed in Sec. |3l the 
distribution function will be exactly modified by the source term 

5f = ---^ t f. (5.1) 

Now we define a corresponding continuous central moment of order (m+n) due to change 
in the distribution function as a result of a force field as 

/CO /-CO 
/ ^ F (^- M ,n£y-"y«^y- (5-2) 
-CO J —CO 

Substituting Eq. fl5.1|) in Eq. (J5.2|) and integrating by parts by making use of the asymptotic 
limit assumptions lim^ ±co (^ x -u x ) m f(x,y r ^ Xr ^y)=0 and lim^ ±0O (^y-Uy) n f(x,y^ x ^y) = 
0, for m,n>0, we get 

p /-CO /-CO 

f£ = m-^/ / f(Z x -U X ) m - 1 (Zy-Uy) n dZ x dZy + 

P J -co J — CO 

P 



P f' CO /-CO 

»-/ / /(&-« X ) m (gy-«y) n " 1 ^x^y. (5-3) 

jO J co, J co 
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From the definition given in Eq. (|4.1|), Eq. (|5.3[) reduces to an exact identity between con- 
tinuous central moment of the source term of a given order to the components of the 
continuous central moment of the distribution function of an order lower acted upon by 
a force field: 

and for the special case of the zeroth central moment of the source as Tq = 0. This is a 
key result based on which the rest of the derivation follows. Thus, we can enumerate the 
exact values of the central moments of sources in an increasing order as 



1 


= o, 


fF 

1 X 


= -n , 
p 


fF 

y 


— — 1*0/ 


fF 

XX 


= 2-flx, 
P 


fF 

yy 


= 2^h y , 
p 


fF 

^y 


= ^n y +^n x , 
p p 


fF 

xxy 


= iEin xy +^n 
p p 


fF 

xyy 


p 33 p 


fF 


F X ft Fy 

= 2-\\ x:i; , ■ 2- 

p 33 p 



(5.5) 



xy, 



'■xxy- 



Note that if we set H x < n y n =n^ v « in Eq. 05.5D , i.e. ignore non-equilibrium effects, we arrive 
at the the derivation given in [63] as a special case. 



6 Cascaded Central Moment Lattice-Boltzmann Method with Forc- 
ing Terms 

Defining a discrete distribution function supported by the discrete particle velocity set lt a 
as f = \f a ) = (/o,/i ,fi, ■ ■ ■ ,fs ) T / an d a cascaded collision operator as O c = | d c a ) = (Cl c , Cl\, Cf 2 , . . . , Qg ) 7 
as well as a source term as S= \S a ) = (So,Si,S2,---,Ss) T based on central moments, we ob- 
tain the lattice Boltzmann equation (LBE) as a discrete version of Eq. (|3.1|l by temporally 
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integrating along particle characteristics as follows Il63ll : 

f^ + ^ a/ t + i)=f^ r t)+n^ >t) +J t t+ \ ailt+ ^ M+d) de. (6.i) 

Here, the collision operator is written as 

n»=ns(f / g) = (x;.g) B/ (6.2) 

where g= \g a ) = (go,gx,g2,---/g8) ■ The hydrodynamic fields are obtained from the distri- 
bution function as 

P = E/« = </«IP>/ (M 

a=0 
8 

P u i = Y^fAi=(fxKi)4&x,y. (6.4) 

For improved accuracy in recovering Navier-Stokes solution, using a semi-implicit rep- 
resentation for the source term, i.e. the last term in the above equation (Eq. (|6.1[l) as 

It t+ls cc(lt+tM+e) de = l [ S <Y(l?,f) + S «(^+^+i)] ' so that E q- ED is made effectively ex- 
plicit by using the transformation f a =fu — ^S a to reduce it to |63| 

7«(^+^*+i)=7«(^/0+n« ( ^)+s a( ^ t) . (6.5) 

It may be noted that, as in II59II , we first represent collision as a cascaded process in which 
the effect of collision on lower order central moments successively influence those at 
higher orders in a cascaded manner. That is, in general, g a =g a (f,gp),j5 = 0,l,2,...,a. — l. 
Furthermore, the form of the source term is derived to rigorously enforce Galilean in- 
variance. The explicit expressions for S a and g will be determined later in Secs.[9]and [lOj. 
respectively. Since the main focus of this work is on improving the collision (including 
forcing) step with features independent of inertial frames, we have only considered the 
standard discretization for the streaming operator. However, as discussed in the Intro- 
duction, other types of discretization schemes could be considered to improve the order 
of accuracy. 



7 Various Discrete Central Moments and Galilean Invariance Match- 
ing Principle 

For determining the structure of the cascaded collision operator g and the source terms 
S tt , we first need to define the following discrete central moments of the distribution func- 
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tion, Maxwellian, and source term, respectively: 



^_jfa.{^«.x U x ) {s tt y Uy) — ({Sax U x ) {t^ccy Uy) \fa}/ 



(7.1) 



a. 




y~!./a^ ( e ccx u x) { e ccy u y) — {{^ax u x) { e ay u y) \fu^)r 



(7.2) 



^^S a (c ax U x ) {c a y Uy) — {{?ax Ux) {Cay Uy) |S a ). 



(7.3) 



a 



where the exact expression for the discrete /^ is not yet known, but can be determined 
as a result of the derivation discussed later. To maintain physical consistency at the dis- 
crete level, we now equate the discrete central moments of the distribution function, the 
Maxwellian and the source terms, defined above, with their corresponding continuous 
central moments, whose forms are known exactly. That is, according to this matching 
principle 



In particular, the discrete central moments of various orders for both the Maxwellian and 
the source terms, respectively, become 




(7.4) 
(7.5) 
(7.6) 




■xxyy 



o, 
o, 

c 2 s p, 
c 2 s p, 

0, 




(7.7) 
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and 





= o, 


<r x 


= —K , 
P 




= —K , 
P 


O-xx 


= 2—K x , 
P 


°yy 


Jy~ 

P J 


&xy 


F x ~ F y ^ 
— Ky~\ K x , 
P P 


vxxy 


rjF x — .Fy — 
P P 


a xyy 


Fx— \^y— 

— Kyy-T^. K X y, 
p JJ p 3 


Vxxyy 


= 2 — K x yy-\-1 — K X 



We also define a discrete central moment in terms of the transformed distribution function 
to facilitate subsequent developments as 

= YJoc ( e Kx-u x ) m (e a y-u y ) n = ( {e ax - u x ) m {e ay - u y ) n \f a ) . (7.9) 

a 

Owing to the transformation discussed in Sec. it follows that 

- - I- 

K x m y n = K x m y n — -(X x m y n . (7.10) 



8 Relation Between Various Discrete Central Moments 

Equation (|7.8|) is given in terms of the discrete moments of the original distribution func- 
tion f a . However, the cascaded central moment LBM with forcing term provides evolu- 
tion in terms of transformed distribution function / (Eq. I|6.5[)). Thus, it is important to 
write all the expressions in terms of the central moments of f a , or, equivalently, K x m y n . 
Thus, by recursive application of Eq. (|7.10|) using Eq. (|7.8|) to successively higher orders, 
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we get the following exact relations up to the third-order central moments as 



K = K 0/ 

- 1 F x ~ 

Kx = K x + - — K , 
2 P 



K V 



~ F x ~ \Fl~ 

K xx — K xx -\ K x -\- — ^tKq, 

P 2p 2 



Fx^ l F y, 
7 Ky + 2^ 



K yy ~ K yy + ~ K y + n~T2 K Q' 



^ 1F X ~ lFy~ lF x Fy^ 

K xy - Kx y + 2~p Ky 2~p Kx 2~p rK °' 

~ F x ~ lF y ~ F x F y ~ lFl~ 3 F^Fy^ 
Kxxy - K xxy+—K X y+^—K xx + -^-K x + --^Ky + --^-Ko, 

Fy^ 1F X ^ F x F y ^ 1F^ 3F X F^ 

K X yy - K xyy + ^K xy + ^ ^K yy + ^ K l/+ 2j0 2 K *+4 pi K °- 



That is, the central moment of the distribution function of a given order can be written 
as a function of the central moment of the transformed distribution function of the same 
order and successively lower orders as well. This can be further simplified by consider- 
ing the three of the lowest order central moments, i.e., conservative moments, which by 
definition are £o = p, k x = K y = 0. This, in turn, leads to Ko = p, k x = — 1 / 2F X , K y = — 1 / 2F y . 
As a result, we have the following relations for the non-conserved central moments up to 
third-order: 



K XX 


— K xx , 


K yy 


= *yy 


K X y 


= K X y, 



K XX y — K 



K X yy — 



F X ~ 1 Fy^ 

xxy~^~ Kxy-I-^ K xx , 
Fy_~ IFx- 



^xyy 



16 



Thus, we can finally write the central moments of the source term in Eq. (|7.8[) in terms of 
the central moments of the transformed distribution function as 



0-Q 


= o, 




= F x , 


Vy 


= Fy> 


"xx 


= o, 


vyy 


= o, 


Cxy 


= o, 




= 2^ 


&xxy 






P 



(8.1) 



Kxyi ^xxi 
V P 

&xyy — ~~^~^yy' p~^ x y' 

a - ?— ? + +5? i a FxF Vjf 

u xxyy — ^xyj/T^ ^ <^xxy\^p_ xx 'pi yy pi X V' 



Thus, higher order non-equilibrium effects in K x m y n and non-linear effect in F£ Fy are ev- 
ident for the central moments of the source terms that are third- and higher orders. Let 
us now explicitly write the central moments of the transformed discrete Maxwellian by 
means of Eq. (I7.10D using Eqs. (|7.7D and d8.ll) to yield 



K 


= P> 




*x 


- _If 

— 2 r x> 




K y 


- _I f 

— 2 y 




~M 
K xx 


= c 2 s p, 




~M 

K yy 


= c 2 s p, 




K xy 


= o, 




^xxy 


p Kxy 


Fy^ 

2p Kxx ' 


%x yy 


2p Kyy 


Fy^ 
p Kxy ' 


K xxyy 


4 Fx 
= c sP -- 


- F y^ 
Kx yy p K - 



(8.2) 



F 2 p2 p c 

— ir / i-ic 

xxy 2p 2 2p 2 yy p 2 y ' 

The main idea in the determination of the collision operator for the cascaded version 
of the central moment method is to relax the central moments of the transformed dis- 
tribution function to its corresponding local attractor, successively at various orders as 
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given in Eq. (|8.2|) (see Sec.[lO]). Before proceeding further to do this, we first need certain 
quantities in the rest or lattice frame of reference, i.e. the raw moments, in which the 
computations are actually performed. These are obtained in the next section. 



9 Various Discrete Raw Moments and Source Terms in Particle 
Velocity Space 

The raw moments, i.e. those in the rest frame of reference, can be related to the central 
moments by means of the binomial theorem [85 . 86 1. For any state variable (p supported 
by the discrete particle velocity set, the transformation relation between the two reference 
frames is thus given by |i63| 



((e« x -u x r(e ay -u y y\cp) = + < 



\<p)+ 



-ay 



£crcr(-i) 1 " 
£crc- i (-i)x 



Y2cfep(-iy4 



W) (9-1) 



where Cq =p\/ (q\(p — q)l). We now define the following notations for depicting various 
discrete raw moments, based on which an operational LBE will be devised later: 



K x myn 


— V f e m e n — 

a 


Ipin e n \f \ 
v-ax'-ay\j a I i 


K x m y n 


= VI e m e" = 

a 


(e m e n 17 ) 

\'-ax'-ay\J al ' 


a x m y n 


— r<; c"/ - 


\ c &x'-ay ■ 



(9.2) 

(9.3) 
(9.4) 



Note that the superscript "prime" (') is used to distinguish the raw moments from the 
central moments that are designated without the primes. Here, analogous to Eq. (17.1011 , 

the relation K x m y n = K x m y n — j^' x m y n holds. Let us first find expressions for K x ,„ y n = (f a \e™ x e". y ) 
to proceed further. As in |63|, for convenience, we define the following operator acting 
on the transformed distribution function / in this regard: 



<fa 1 +fa 3 +fa 3 + ---)+Hf fil +f fi2 +ff i3 + ---) + --- = 
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where A = {ot-iAiAzr ■ •}/ B = {^1,62/^3,- ■ ■},■■■■ For conserved basis vectors, we write 
them in terms of the hydrodynamic variables and force fields as 

% = (fJp) = Lf«=P, (9-6) 

^ 8 - 1 

k * = =E/« e «*=P M *-9 F *' ( 9 - 7 ) 

a=0 Z 



1 



</>*y> = E f^y=P u y ~ ( 9 - 8 ) 

a=0 

and, for the non-conserved basis vectors, using Eq. 09.5b in terms of subsets of particle 
velocity directions as 

8 / A 3 \ 

(9.9) 
(9.10) 





(7. 


k 2 ) = 

raj/ 


8 /A 3 \ 

E7a4= E 

a=0 \ ci J 


*yy = 


(7. 


Fay/ 


8 _ 2 M\ 

E7a e ay — I E ) 
a=0 \ a / 


K X1/ = 


(7 a 


l^ax^cey, 


8 _ 

> = E/a eaxea y = 

a=0 


^xxy 


(7 a 


l e ax e ay, 


8 _ 

^ = E/a e <*x e ay = 
a=0 


K xyy ~ 


(7 a 


e ax e ay, 


8 

) = E/a eaxe «y = 

a=0 


K xxyy = 


(7 a 


1 2 2 ' 
l^ax^ay, 


8 

' 7 a'-ax'-ay 
a=0 



As B s 



a a 
'A 6 B 6 ' 



E-E )®f«> ( 9 - 12 ) 

a 1 y 

E-E ^ ( 9 - 13 ) 



a a 
'A 8 \ 



E ®/* (9-14) 



where 



A 3 = {1,3,5,6,7,8}, 
A 4 = {2,4,5,6,7,8}, 
A 5 = {5,7},B 5 = {6,8}, 
A 6 = {5,6},B 6 = {7,8}, 
A 7 = {5,8},B 7 = {6,7}, 
A 8 = {5,6,7,8}. 

Now, we transform the central moments of the source terms (Eq. (17.81) ) to the cor- 
responding raw moments by considering Eq. (|7.3|) and using the frame transformation 
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relation (Eq. (|9.1jl ). This yields 



Vo = {S a \p}=0, 


(9.15) 


0~ x ~ {S« \&ax) = Fx/ 


(9.16) 


^y = W« 1 *-«y ) = Fy/ 


(9.17) 


^xx = (^a e ax) = 2Fx-Uj, 


(9.18) 


^yy = (^al^ay) = 2F x Uy, 


(9.19) 


^xy = W« \^ax^ay) = FjMy + FyU x , 


(9.20) 



^xxy — (^« l e <uAy) — 2fx I — 



v.vy 



-Mx"y 



^xyy — (S« | e ax^a V ) — Fx 



2 

ay/ 



vyy 



+ "y +2Fy 



^xy 



+ U x Uy 



(9.21) 
(9.22) 



'xxyy 



(S«|<&<£ y >=2F*uJ 



-2F yUy 



y + " 2 x 



2F X ~ 2F V ^ 



F 2 , 



F y- 

p **y^^2 K yy + ^2 K **~ 



F* M y ^y 



(9.23) 



Clearly, the raw moments of source terms for third-order or higher contain non-equilibrium 
and non-linear contributions. Eqs. (19.21l) - (19.23l) require explicit expressions for central 
moments of transformed distributions such as K xx , K yy , K xy , K xxy and K xyy , in terms of raw 
moments to facilitate computation. They can be readily obtained in terms of raw mo- 
ments from their respective definitions and by using the binomial theorem (Eq. (|9.1ll ) and 
subsequent simplification as 



K xx 



K 



xy 



K xx + F x u x -pu 2 x/ 

— i r 2 
K yy< t y u y~P u y 

J 1 

K xy + -{F x u y + F y u x )-pu x u y , 



(9.24) 
(9.25) 

(9.26) 



for second-order and 



^yy 



\-xxy 



K xyy 2-U y K X y- 



- UrflL, - -F x u y -F y u x u y +2pu x u 2 y , 



yy 



K xxy 2-U x K xy U y K xx ^ V U x 



2 F x U x Uy + 2pU 2 x Uy, 



(9.27) 
(9.28) 



for third-order moments. Based on the above, we now obtain the source terms projected 
to the orthogonal moment basis vectors, i.e. (Kp\S a ), ft = 0,1,2,. ..,8, which would then 
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provide corresponding explicit expressions in terms of the particle velocity space. Thus, 
from Eqs. (I2.5D and (|9.15I) - (I9.23I) , the following projected source moments are derived: 



m s Q = 


(K \S K ) 


= o, 


m\ = 


(Ki\S a ) 


= fx, 


m\ = 


(K 2 \S K ) 


= h> 


m\ = 


(K 3 \S a ) 


= 6{F X 


m% = 


(K*\S K ) 


= 2{F X 


m% = 


(K 5 \S a ) 


= (F x u 


m% = 


(K 6 \S a ) 


= -6F. 



m s 7 ={K 7 \S a ) 
m s s = (K 8 \S a ) 



Ix + FyUy), 



r y u yJ 



-F y u x ) 



+ U x Uy 



-3F V ^ 



-3R ^ 



K X x 
P 



-6Fj^ 



18FrMi 



M/l/ 



+u 



2 |)+18Fy«y 



XX 

p 



F F F 2 F 2 

18—K xy y + 18^K xxv +9^K vy +94rK 



xxy 



36 



F x U y ^ FyUx ^ F x Fy 



n 2 '^yy 



(9.29) 
(9.30) 
(9.31) 
(9.32) 
(9.33) 
(9.34) 

(9.35) 
(9.36) 



(9.37) 



In Eqs. d9.35D - l|9.37 |), k XX/ K yy , K xy , K xyy and K xxy can be obtained from Eqs. d9.24b - J9.28ft , 
respectively. This can be written in matrix form as 



/C T S = (/C-S) a = ({K \S a ),(K 1 \S a ),(K 2 \S a ),...,(K 8 \S a )) 
= {mlM\,m\,...,fhl) T =\rh s cl ). 



(9.38) 



Now, by exploiting the orthogonal property of /C [63|, i.e. /C 1 = /C T -D 1 , where the diag- 
onal matrix is D=diag((K |^o),(Xi|K 1 ),(K 2 |K 2 ),...,(X8|X 8 ))=diag(9,6,6,36,4,4,12,12,36), 
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we exactly invert Eq. (|9.38|) to finally obtain source terms in velocity space S a as 



So = 


1 

9 H&S 


(9.39) 


c 


1 

— [6in\ — fn^ + 9/w| + 6m 7 — 2mg] , 




s 2 = 


1 

— [6m s 2 -m\- 9m% + 6m s 6 - 2m s 8 ] , 


(9.41) 


s 3 = 


1 

— [-6m\-ml+9rh\-6m s 7 -2m\], 


(9.42) 


s 4 = 


1 

— [-6m s z -fh s 3 -9mi-6fh s 6 -2fn s 8 ], 


(9.43) 


s 5 = 


1 

— [6m\ + 6m| + 2m% + 9m| - 3m| - 3m^ + m|] , 


(9.44) 


s 6 = 


1 

— [-6mf + 6m| +2m| - 9m s 5 -3m 8 6 +3m s 7 + m|] , 


(9.45) 


s 7 = 


1 

— [-6mf-6m|+2m|+9m|+3m|+3m^+m|] / 


(9.46) 


s 8 = 


1 

— [6mf - 6m| +2m| - 9m| + 3m| - 3my + m s 8 ] . 


(9.47) 



This is the explicit set of expressions for the source terms in velocity space S a given in 
terms of and K x myn, with 2<(m + n)<3 and < m,n < 2. 

Again, using the orthogonal property of /C, we can obtain the raw moments of the 
collision kernel 



E(K'S)"«, = E<^l« y >^ (9-48) 



which is of central importance in the subsequent derivation. Note that for collision in- 
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variants, g =g 1 =g 2 = 0. We get 

E(^ - g)« e «* = E( K |sl e «)S/8 = 0, 
E( /c -8)* e «y = E( JC i8|e«y)g' J 8 = 0, 

ft 

E( /C -g)-4 = E(^l^)^ = 6^3-2^4, (9.49) 

a p 

^(/C-g) a e ax e ay = E(K/3K z e ay )g / 3 = 4g 5/ 
E( /c, g)« e L e «y = E(^l e ^ e «y)^6 = -4g6/ 
E( /C -g)« e ^ e «t/ = E( K /5l e ^ e «y)^ = _4 #7' 

a /3 

E( /C -g)« e ^i / = E( K /5l e ^ e «v)^ = 8 ?3+4^ 8 - 

Finally, the LBE in Eq. 06.5b can be rewritten in terms of collision and streaming steps, 
respectively, as 

7,(^0 = /^O+n^+s^), (9.50) 

7,(^ + "^+l) = Jj3,t), (9-51) 

where the symbol "tilde" (~) in the first equation refers to the post-collision state. In 
terms of the transformed distribution, the hydrodynamic fields can be computed by 
means of the following: 

P = E7 a =<7», (9-52) 

(Y = 

8 - 1 - 1 

P U i = Ylfc t e «i+i F i=(fc t \ e «i) + i F i' ieX >V- (9-53) 

10 Structure of the Collision Operator: Cascaded Central Mo- 
ments 



Let us now arrive at the expressions for the cascaded formulation of the collision op- 
erator using central moments in the presence of forcing terms based on the results ob- 
tained in the last few sections. The basic procedure can be stated as follows. Beginning 
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from the lowest order central moments that are non-collisional invariants (i.e. k X x and 
higher), they are successively set equal to their local attractors based on the transformed 
Maxwellians (Eq. d8.2[l ). This step provides tentative expressions for g a based on the 
equilibrium assumption. We then modify them to allow for relaxation during collision 
by multiplying them with corresponding relaxation parameters [59 J. In this step, given 
the cascaded nature of the collision (i.e. g a =g a (f,gp),fS=0,l,2,...,x—l, or the dependence 
of higher order terms on those that are lower orders), care needs to be exercised to mul- 
tiply the relaxation parameters only with those terms that are not yet in post-collision 
states (i.e. terms not involving gp,f> = 0,1,2,. — 1 for g a ). Various details involved in 
this procedure are given in [63]. For brevity, here we summarize the final results which 
are as follows: 

h = ^{lp+p( u x + u y)-(%x+% y )-\(v xx +a yy )y (10.1) 

CO4. [ ~f 1 / / "I 

= ^y(ul-u 2 y )-{K xx -K yy )--(a xx -a yy ) \, (10.2) 

g6 = -^<^x u y + K xxy- 2u xK X y-UyK xx --Crxxy>--Uy{3g 3 +gi) 

-2uxg 5 , (10.4) 

C0 7 ( 2 1_ ] 1 ^ , 

g7 = -^<^xUy + K xyy -2 Uy K xy -U x K yy --(7xyy>--U x (3g 3 -g 4 ) 

-2u y g 5 , (10.5) 
o; 8 f 1 „ 2 7 |W n -' ~ -' ?-' ?-' 

88 = t\9 p p x v~r xx vy xKx vv yKxx v + xK yy + y Kxx 

+4u x Uy% y - ~v xxyy \ -2g 3 - ~mJ(3S&+&) - \u 2 x(3g3-g4) 
-4u x Uyg 5 -2uyg(,-2uxg 7 . (10.6) 



t ?4 



In the above, the raw moments of various orders, i.e. K x myn for different m and n are 
required, which may be obtained from Eqs. d9.6D - d9.14ll . Similarly, the raw moments of 
sources of various orders, i.e. d~ x myn needed in the above are given in Eqs. d9.15[) - (|9.23j) 
(see Sec. [9]). Here, cor, where f> = 3, 4,5,..., 8, are the relaxation parameters, satisfying 
< cop < 2. When a multiscale Chapman-Enskog expansion |64| is applied to this cen- 
tral moment LBM based on central moments, it recovers the Navier-Stokes equation with 
the relaxation parameters CO3 = oo x and 0^4 = 0^5 = to v controlling bulk and shear viscosi- 
ties, respectively (e.g., v = cl — |) ) II63II . The rest of the parameters can be adjusted to 
improve numerical stability. 
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11 Cascaded Central Moment Lattice Boltzmann Equation 

The collision step with the addition of forcing terms (see Eq. (|6.2[) and Eq. (|9.50[l ) in the 
stream-and-collide procedure of the LBM, obtained by matching those of the continuous 
Boltzmann equation as discussed in the previous sections, is expanded element-wise and 
can be summarized as follows: 



7o 


= 7o + [£o-4(g3-g 8 )]+S , 








= fl + \go+gl-g3+g4 + 2(g 7 -l 


js)]+Si, 




% 


= f 2 +\£o+g2-g3-g4:+2(g6-!: 


rs)]+s 2/ 




% 


= ?3+[g0-gl-g3+g4-2{g7+^ 


?8)}+S 3 , 




% 


= 7i+\g0-g2-g3-g4:-2{g6+!: 


rs)]+s 4/ 




% 


= 7 5 + lSo+gi+g2 + 2g 3 +g 5 -g ( 


>-g7+g8 


}+s 5 , 


% 


= f6+lg0-gl+g2 + 2g 3 -g5-gt 


i+g7+g8 


]+s 6 , 


% 


= 7 7 +lgo-gi-g2+2g 3 +g 5 +g ( 


i+g7+g8 


}+s 7 , 


% 


= 7 8 +[g0+gl-g2 + 2g 3 -g 5 +g ( 


i-g7+g8 


]+s 8 . 



The collision kernel ge needed here can be computed from the expressions given in the 
previous section (see Sec. [TOjl. The source terms in the velocity space can be obtained 
from Eqs. I|9.39|l - (|9.47|l (see Sec.©. The remaining streaming step is carried out as usual 

by using the post collision values / obtained from above. Once the local distribution 
function is known, macroscopic fluid density and velocity fields satisfying the Galilean 
invariant Navier-Stokes equations in the presence of force fields can be obtained from 
Eqs. (|9.52|) and (|9.53[) , respectively. 



12 Cascaded Collision Operator with Reduced Compressibility 
Effects 

While a main goal of this work is the introduction of a self-consistent approach based on 
the continuous Boltzmann equation to incorporate non-equilibrium effects into the cen- 
tral moment approach for general applicability, it is also useful to consider its limiting 
cases. For example, the incompressible limit of fluid flow corresponds to considering 
very small deviations from the local equilibrium, a special case with various applica- 
tions. In particular, this would allow simple representation of incompressible turbulence 
considered later in this work. 

Being a kinetic approach, the lattice Boltzmann method is inherently compressible in 
nature. On the other hand, when it is desired to reproduce the "incompressible" Navier- 
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Stokes equations as mentioned above, it is important to reduce such compressibility ef- 
fects. An approach in this regard was introduced earlier in (87). Here, we will extend 
it further in the context of the central moment LBM in the presence of forcing. It may 
be noted that the fundamental expressions for the continuous central moments for the 
local equilibrium as well as the forcing given in Sees. H] and |5l respectively, from which 
their discrete counterparts are derived, remains unchanged for this case. However, the 
key element to incorporate a systematic reduction of compressibility effects lies in the 
following careful definition of the raw moments of the hydrodynamic fields: 



P = E/«=E/«' 

x=0 ce=0 



po u 



8 8 i 

Y^fu e <x = Y2fu e a + ^ ~F ^ ■ 

a=0 



(12.1) 
(12.2) 



a=0 



where p = po+Sp. Here po and Sp are the constant reference value and fluctuations of 
density, respectively. That is, in the above, contributions of density fluctuations are elim- 
inated from first-order moments representing the components of momentum. Thus, we 
get 



(7» = P, (12-3) 

(7>«*> = PoUx + 2 F *' < 12 - 4 ) 
- 1 

(/J e «y) = Pou y + -Fy. (12.5) 



Using the procedure discussed in the previous sections and with the above specialized re- 
definition of the conserved moments, we obtain, after some simplification, the cascaded 
collision operator with reduced compressibility effects. They are reported here in the 
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following: 



g3 
gi 

g5 
g6 

g7 



0^(2 

12 1 3 

Cl'4 

T 



^-\(po-Sp)u x Uy-K 1 



p+(p -5p)(u 2 x + u 2 J )-(K xx + K yy )--(a xx +a yy ) }, 
(po-Sp) (u 2 x -u 1 )- (k xx -k ) - - (a xx -a ) \ , 



x y 2 axy \' 



X \ ( 2 P0-^p)u 2 x U y + K xxy -2u x K xy -U y K xx --(T, xl , 



-2u x g 5 , 



007 



— \ (2po - Sp) u x ii; + K xyy - 2u y K xy -u xKyy - -a : 



yy 2 



'xyy 



'2 U y( 3 ?3+?4) 



-2 U ^g3-g^) 



-2u y g 5 , 

^{lpH3poSp)uW 



(12.6) 
(12.7) 
(12.8) 

(12.9) 
(12.10) 



y 



wyy 2-U x K xyy 2UyK xxy + U 2 K yy + U y K xx 



+Au x U y K xy 



1^ 

' 2 axx yy 



1 1 

fa - -^u 2 y (3g 3 +gi) - -u 2 x (3g 3 -gi) 



Au x u y g 5 -2ii y g 6 -2u x g 7 . 



(12.11) 



The above collision operator that selectively introduces density fluctuations where nec- 
essary can reduce compressibility effects for a inertial frame invariant flow field while 
retaining its linear acoustics. It can thus allow for a better comparison with "incompress- 
ible" macroscopic fluid dynamic equations, particularly for turbulent flows as discussed 
later. 



13 Factorized Central Moment Model for Collision 

In this section, we will derive an alternative form of the central moment LBE with forcing 
terms based on a different choice of the local attractor in the collision step for improved 
physical coherence. Continuous kinetic theory, as originally initiated by Maxwell |88|, 
features two important properties for the local equilibrium or the Maxwell distribution 
- Galilean invariance and factorization in Cartesian components of the particle velocity. 
As discussed recently fl39U62L it could prove useful to inherent these properties at the dis- 
crete particle velocity level. The use of central moments maintains Galilean invariance 
by construction. Factorization property of the distribution function implies that particle 
velocities are random variables. An extension of the factorization idea beyond equilib- 
rium as a model for describing the relaxation process during collision was proposed to 
construct local attractors II62II . Specifically, the basic postulate behind this model is that 
the Cartesian products of the post-collision values of the orthogonal central moments of 
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lower orders that are not in equilibrium forms as the basis for the attractors of the higher 
order moments. Here, we further extend this to include source terms so that the model 
can incorporate force fields. Thus, the attractors for central moments of different orders 
are given as 



^at 

K y 

<^xy 

^xxy 

K at 
^xyy 

K at 



0, 

o, 



K xKyy 
KxxKyy, 



o, 

= 0, 
= 0, 



(13.1) 
(13.2) 

(13.3) 
(13.4) 
(13.5) 
(13.6) 



while the second-order longitudinal central moments are obtained from Maxwellian as 
given in an earlier section (see Sec. 0, i.e. k xx = K%t = pc 2 . In essence, the distinguishing 
feature of the factorized central moment lies in the use of modified attractors for third 
and higher order moments. Now, using the following central moment identity of the 



post-collision state K x m y" — K x 
m = 0,n = 2, we get 



-E/3(^l(e a x-Wx) m (e a y-%)")^, for m = 2,n = and 



Kxx = K« + (6g 3 +2g4)/ 

Kyy = Kyy+(6g 3 -2g 4 ). 



(13.7) 
(13.8) 



Note that it also follows that K xx = k xx and Kyy = Kyy. We can then rewrite everything in 

terms of transformed raw moments, i.e. k xx = k xx — pu x +F x u x and Kyy = Kyy— ptiy+ FyUy. 
These yield 



v yy 



K xx - P u l +FxU x + {6g3 +2g4), 

/ 

^yy ~pU 2 y + FyUy + (6g 3 ~ 2g 4 ) • 



(13.9) 

(13.10) 



In effect, the attractor for the fourth-order moment, i.e. Eq. (113.61) reduces to 



xxyy 



K v 



-pu 2 x + F x u x + (6g 3 +2g A ) 



Kyy ~ P u y + F y u y + (6g3 ~ 2#4 ) 



(13.11) 



Now, to obtain an operational step in terms of the transformed variables, we use the 
2&x m y n to finally get the following expression for the fourth-order 



relation K x m y n 
central moment 



.exit 

■ K- X my,i 



K 



xxyy 



K XX ~PU 2 x + FxU x + (6g 3 + 2f 4 ) J \Kyy ~ pu] + FyUy + (6g 3 - 2g 4 ) 
F _ F... F 2 . _ F2 _ F..F.._ 
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By replacing K y m,,„ 



c x m y" 



with K x myn as given above, we can derive the collision kernel with 
factorized central moments as attractors. It follows that the expressions in Sec.[lO]for ga, 
ft = 3,4,5,6,7 are the same as before with the exception for g$. The expression for g$ in 

Eq. 010.6)1 is modified such that term \p{=K xxyy ) in this equation is now replaced by ^ xxy y 
given in Eq. (113.12b . In a similar vein, the above expression can be modified for reduced 
compressibility effects (see Sec. [12]) as 



xxyy 



*xx ~ (po ~ Sp)u 2 x +F x u x +(6g 3 +2gi) 

Kyy ~ (PO ~ &P) Uy + F yUy + (6g 3 ~ 2g 4 ) 



Fx 

K 



xyy' 



F F 2 

ry^ y ^ 



F 2 ~ 

2p i K yy L 



FxFy 



2 Kx v 



(13.13) 



to modify gs in Eq. dl2.HI) . 



14 Galilean Invariant Filtered Turbulence Representation using 
Lattice Kinetic Framework 

Based on the various elements derived in the previous sections, we are now in a posi- 
tion to construct an approach for simulation of Galilean invariant turbulent flow field by 
incorporating appropriate turbulence models in the LBM. The starting point in the statis- 
tical continuum description of turbulence is the Reynolds decomposition of the velocity 
field of the fluid into 'resolved' and 'unresolved' parts. The resolved part is obtained 
by applying either some averaging in space or time (in the Reynolds Averaged Navier- 
Stokes (RANS) context) or by applying a filter (in the LES). Application of this decompo- 
sition to the Navier-Stokes (NS) equation leads to additional unknown terms involving 
products of the unresolved fields, which are Reynolds stresses (in RANS) or the subgrid 
stresses (in LES). This closure problem then becomes the main focus of turbulence mod- 
eling. Due to the scale invariance property of the NS equations II83II , the averaged and the 
filtered equations, as well as the additional stress-like closure terms have similar forms. 
Thus, a unified statistical approach may be adopted for turbulence modeling. It is inter- 
esting to note that ideas based on kinetic theory provided the original inspiration for the 
Reynolds decomposition [89| as well as early works on developing turbulence models. 

The underlying motivation here is to develop a unified statistical averaged descrip- 
tion (for RANS) or formal spatial filtered representation (for LES) of inertial frame invari- 
ant turbulence in a kinetic approach based on the LBM derived in earlier sections. This 
would also allow reconciliation of continuum and non-continuum effects on turbulence 
as discussed recently [77J. The following notation for Reynolds decomposition is adopted 
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here. For any scalar <p, vector v and tensor TW, we have 

<p = <p+<p , with tp = 0, 



v + u , with p =0, 
Tij+T^ with T-=0, 



where ( • ) is an operator representing either some form of statistical average or filter to 
obtain the resolved part and the symbols with primes denote the unresolved parts of the 
field. As discussed in f77\, application of the above decomposition directly to the con- 
tinuous Boltzmann equation (Eq. d3.ll) ) leads to certain difficulties. In particular, using 
f = f+f for the distribution function in Eq. (|3.1)) , which leads to a statistically aver- 
aged kinetic equation, does not provide a clear interpretation of turbulence physics. The 
local collision term needs to model all essential physics, including the non-linear and 
non-local momentum transfer effects of turbulence. Moreover, the use of the averaged 
attractor based on the Maxwellian f M within the collision term Q (/,/) leads to model- 



ing difficulties since exp 



~2cf 



^ exp 



. Thus, an alternative approach is 



needed to coherently represent continuum /kinetic effects on turbulence. 

To circumvent these issues, a transformation for the velocities was recently suggested f77\ , 
which is adopted here to provide Galilean invariant turbulence representation. The key 
element is to clearly separate the advective turbulence effects due to unresolved veloc- 
ity field u from the dissipative collision that represent microscopic effects. This is ac- 
complished by an inspection of the local Maxwellian given in terms of the microscopic 
particle velocity g and the macroscopic fluid velocity it. That is, it consists of the term 
involving the peculiar velocity 

as its argument which should be made independent of the unresolved part of the macro- 
scopic fluid velocity u, when the averaging operator is applied. That is, 



■It 



should be transformed appropriately, which can be accomplished by defining a new vari- 
able ~rf as 

■f = ~t-lt'. (14.1) 
Now, the Maxwellian in the transformed peculiar velocity if — JL commutes with the 



operator for averaging or filtering. That is, exp 



2Z| 



exp 



. This facil- 



itates the separation of various aspects of turbulence physics modeling. Based on this a 
new distribution function h(x,lf,t) and its local Maxwellian are defined by 



htf,lt,t)=f&,~t t t), 



(14.2) 
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(14.3) 



and _^ _^ 

respectively. The continuous Boltzmann equation, i.e. Eq. d3.ll) (without the forcing term 
for simplicity) is then transformed into a modified kinetic equation in terms of rj and h 
as follows. From Eq. (|14.2|l , ^ v h=^zf. When rj = constant, we have x ^) v = ^ ' x lt' 
and (d t t) n =dt~u . Hence, the derivatives in new variables are 

d t h = d t f+(d t ~t) r ^zf = dtf+(drf') r ^ n h. 
The continuous Boltzmann equation is thus modified to [77] 

d t h+jf- ^ x h+lt' ■ ^ x h - ~ct' ■ ^ v h = Q(h,h), (14.4) 

where ~ct' = dt it' -\-tj- V x lt ' +lt' ■ V x lt' . Considering incompressible flows, where the 
unresolved velocity field satisfies V • ^'=0, we get it' ■ v x h=v x (hlt') and it '■^ x h = 
^ x (hlt'). As a result, Eq. (|14.4[) is further simplified to 

d t h+lf x h+\t x (hlt')-^ x {h~a r ') = Ci(h,h). (14.5) 

Now, applying the statistical averaging or filtering operator on Eq. (114.51) and using 
d(h,h) =Cl(h,h), we get the new kinetic equation 

d t h + if ■ ^ x h + ^ x ( hit' ) - ^ x { hit' ) = Q (h,h) . (14.6) 

The averaged density and momentum can then be obtained by taking moments of h. 
That is, p = Jhdlf ,pu ' — Jhlfdlf . Now, h{t]) using Eq. (|14.6[) is better suited to represent 

turbulence physics for the following reasons. The term involving hit' represents trans- 
port in physical space, i.e. redistributes h to smoothen any gradients, hit' represents 
transport in velocity phase space, and acts as a source /sink for energy cascade [77]. In 
particular, both these quantities can be directly related to continuum based closure mod- 
els. On the other hand, the role of collision operator is then to simply represent averaged 
effect of irreversible molecular collisions. The averaged kinetic equation, Eq. (|14.6|) , can 
be further simplified by considering the following simple microscopic closure II 7711 

hit' « hlt' = 0, (14.7) 
h~tt' « h^' =h^ Y -( lt'lt' ) (14.8) 

That is, h is uncorrelated with both u ■ and a-, which reproduces the averaged momentum 
equations with additional Reynolds stress terms u-u- that can be closed by means of any 
conventional macroscopic turbulence models. Thus, Eq. (114.6b can now be rewritten as 

d t h+ ~f ■ ^ x h= fi(M) + ■ ( it' it' ) ■ ^ v h, (14.9) 



31 



which represents the evolution of the statistical averaged /filtered turbulence field by 
means of the Reynolds stresses that appear as a forcing term in a kinetic framework. 

Let us now develop a Galilean invariant lattice kinetic equation, i.e. which provides 
inertial frame invariant representation with respect to the resolved velocity field obtained 
by statistical averaging /filtering. For brevity and to avoid the use of additional new 
notations, let us rewrite Eq. (|14.9|) by replacing h by / (and Jj and £) to make use of the 
developments of the previous sections. That is, 

d t f+ 1 ■ ^*/=0(/,/) + ■ ( it' It' ) • (14.10) 
from which the resolved hydrodynamic fields can be defined as follows: 

p = J fit. Pjt_= J f?df. (14.11) 
Now, we define operator averaged continuous central moments as 

/CO rCQ 
/ ftfx-Ux) m (Sy-Uy) n dS x d$y (14.12) 

~ M 

and similarly for the continuous central moments of the Maxwellian Tl x m y n based on re- 
placing h M and if by f M and g, respectively. The Cartesian components of the unre- 
solved turbulent Reynolds stresses may be written as 

4 = -d x ( u x u x ) -dy(u' x u' y ), (14.13) 
^ = -d x (u x u y )-d y (u y u y ), (14.14) 

where = (a' x ,a'), from which we analogously define a source/ sink continuous central 
moment as 

/OO /»0O / 
/ Sf (k-ttxHSy-Ky)"^^. (14.15) 
-ooj — oo 

Here, 5f a = —_a_-V^f. It readily follows from Eq. d5.4|) that T_ x m y n also satisfies the 

following exact identity T^ x m y n = ma' x Tl x m-i y ii+na' y Yl x m y n-i. That is, the statistical aver- 
aged/ filtered central moment of sources/ sinks due to unresolved fields of a given order 
is dependent on the product of the Cartesian components of the gradients of turbulent 
stresses with the lower order central moments of the averaged/ filtered distribution func- 
tion. The corresponding discrete central moment LBM can be devised by considering the 
following averaged representation of discrete vectors supported by the particle velocity 
set: f = = {fj d J_ v f_ 2 >---j_j & ) T , g= If,) = %4.i'h'---'h )T ' - = l^s) = (SoA,S2,...,Sg) T , 
and n c = O c (f,g) = (/C-g) = (n^Q^Q^,. . .,Og) T , and invoking Galilean invariance match- 
ing principle, i.e. matching the continuous and discrete central moments of various quan- 
tities at successively higher orders as discussed in earlier sections. In particular, the sta- 
tistical averaged /filtered discrete collision operator O c can be obtained by considering 



32 



reduced compressibility effects and factorized attractors as in Sec. [131 Furthermore, the 
corresponding source terms in velocity space S can be constructed using the procedure 
outlined in Sec. [9] The operator averaged LBE, in terms of the transformed distribution 
function / for improved accuracy, can be finally written as 

J^(lt + ^ K 4+l)=l a (^A)+Cy a{ltA) +S a{1 t, t y (14.16) 

where f_ a =f a ~ \S X . Here, as before, we have adopted the standard discretization for the 
streaming step (see the comment following Eq. 06.5D ). The resolved hydrodynamic fields 
in the reduced compressibility formulation can then be obtained as 

8 _ 1 , _ 1 

poUi = EZ«^+9£i = ^J ea ') + ?£i- iex 'y (14 - 18) 

This provides a minimal lattice kinetic equation for incorporating turbulence models, 
where the unresolved turbulent motion are inertial frame invariant with respect to the 
resolved fluid motion. Here, we clarify the meaning of this statement as follows. Unlike 
other areas in fluid mechanics, where models have been developed starting from contin- 
uous kinetic theory, its role for fluid turbulence has been more limited. This is mainly 
due to the fact that kinetic theory generally considers distinct scale separation of physical 
processes. On the other hand, turbulence is a flow phenomenon intrinsically containing 
a continuous spectrum of scales with no scale separation. As such, therefore, turbu- 
lence modeling developments have to rely much on phenomenology whose mathemati- 
cal forms are then constrained by invariance principles (e.g. material frame indifference 
and inertial frame invariance mentioned earlier in the introduction) and realizability con- 
siderations 1 81119011 . Thus, except for some early models such as those based on mixing 
length concepts and derivation of some recent phenomenological models (e.g. [54j) based 
on kinetic theory, turbulence modeling developments are generally based on macroscopic 
models. The ultimate goal of our central moment approach for the filtered kinetic equa- 
tion discussed above, is, then, to simulate resolved turbulent fields which are inertial 
frame independent, when an appropriate macroscopic turbulence model for the unre- 
solve field is used in the forcing term. 

15 Summary and Conclusions 

A discrete lattice kinetic model for the continuous Boltzmann equation, including forc- 
ing, based on central moments is derived. The collision operator as well as the source 
term of this lattice Boltzmann equation is constructed by matching the corresponding 
continuous and discrete central moments successively at various orders. The local at- 
tractor of the collision operator is constructed to satisfy the factorization property of the 
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Maxwellian during relaxation process. An exact hierarchical identity of the central mo- 
ment of sources, that incorporates non-equilibrium effects, is maintained at the discrete 
level. The resulting approach provides Galilean invariant hydrodynamic fields in the 
presence of any external or self-consistent internal forces in a discrete kinetic framework. 
It is further extended to incorporate reduced compressibility effects for better represen- 
tation of incompressible flow, a limiting case. An important physical characteristic of 
turbulent flows is that it is inertial frame independent for all or any subset of scales. In 
consequence, for general applicability, all turbulence models and their simulation ap- 
proaches, should satisfy this requirement. A statistical averaged / filtered lattice kinetic 
equation based on central moments that maintains Galilean invariant representation of 
unresolved fluid motion with respect to the resolved fields of turbulent flow is thus con- 
structed. The formalism presented here can extended to other lattice velocity sets and in 
three-dimensions as well as to other physical problems such as complex fluids. 

In this regard, we make the following remark on the development of more efficient 
schemes for the former aspect. Symmetry and finiteness of the standard lattice sets lead 
to degeneracies of higher order moments in terms those at lower orders that can result 
in frame dependent contributions to viscous stresses. This necessitates considerations of 
extended lattice sets. In this case, it is proposed that the central moments relaxation (as 
well as forcing) be considered only up to those higher order moments that have bearing 
on the physics of hydrodynamics, such as stress tensors and heat flux vectors. In turn, this 
would impose Galilean invariance of the macroscopic description of the fluid motion. 
The relaxation of the rest of the higher moments (including forcing) related to the fast 
kinetic or ghost modes can be considered in terms of the standard or raw moments. Here, 
the form of the hierarchical identity for the sources derived in this paper for those higher 
(kinetic) moments would be the same with the simple replacement of the central moment 
terms by the corresponding raw moments. It is envisaged that such mixed central/ raw 
moment approach would exhibit interesting mathematical structures for the resulting 
collision operator and the sources, as well as being computationally more effective. This 
strategy is currently under investigation. 
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